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LES on unstructured deforming meshes: 
towards reciprocating IC engines 

By D. C. Haworth 1 AND K. Jansen 2 


A variable explicit/implicit characteristics-based advection scheme that is second- 
order accurate in space and time has been developed recently for unstructured 
deforming meshes (O’Rourke Sz Sahota 1996a). To explore the suitability of this 
methodology for large-eddy simulation (LES), three subgrid-scale turbulence mod- 
els have been implemented in the CHAD CFD code (O’Rourke & Sahota 1996b): 
a constant-coefficient Smagorinsky model, a dynamic Smagorinsky model for flows 
having one or more directions of statistical homogeneity, and a Lagrangian dynamic 
Smagorinsky model for flows having no spatial or temporal homogeneity (Meneveau 
zt al. 1996). Computations have been made for three canonical flows, progressing 
towards the intended application of in-cylinder flow in a reciprocating engine. Grid 
sizes were selected to be comparable to the coarsest meshes used in earlier spectral 
LES studies. Quantitative results are reported for decaying homogeneous isotropic 
turbulence, for linear (non-solenoidal) strain of homogeneous isotropic turbulence, 
and for a planar channel flow. Computations are compared to experimental mea- 
surements, to direct-numerical simulation (DNS) data, and to rapid-distortion the- 
ory (RDT ) where appropriate. Generally satisfactory evolution of first and second 
moments is found on these coarse meshes; deviations are attributed to insufficient 
mesh resolution. Issues include mesh resolution and computational requirements 
for a specified level of accuracy, analytic characterization of the filtering implied by 
the numerical method, wall treatment, and inflow boundary conditions. To resolve 
these issues, finer-mesh simulations and computations of a simplified axisymmetric 
reciprocating piston-cylinder assembly are in progress. 


1. Introduction 

Contemporary three-dimensional time-dependent models for flow and combustion 
in reciprocating IC engines are based on solutions to Reynolds- averaged governing 
equations (‘RANS’ based modeling; Amsden et al. 1989, Haworth et al. 1990). In 
RANS, the local instantaneous value of a computed dependent variable represents 
an ensemble- or phase-average over many engine cycles at a specified spatial loca- 
tion and crank phasing. In general, two-equation ( k — e) closures have been used 
to model turbulent transport, with standard equilibrium wall functions. Shortcom- 
ings of RANS models have been documented by several generations of turbulence 
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researchers. Discussion of engine-specific issues can lx* found in the reviews by El 
Tahry Sc Haworth (1992, 1996). 

An alternative to RANS is large-eddy simulation (LES). Here the governing equa- 
tions are spatially filtered rather than ensemble averaged. Explicit account is taken 
of flow structures larger than the filter width, which is on the order of the mesh 
spacing, while the influence of unresolved scales is modeled using a subgrid-scale 
(SGS) model. Because statistics of small-scale turbulence are expected to be more 
universal than those of the large scales, LES offers the promise of wider generality 
and reduced modeling uncertainty. 

LES is particularly appealing for IC engine applications. Turbulence model for- 
mulation and calibration traditionally have been carried out in statistically station- 
ary and/or homogeneous flows for simple geometric configurations. To bring these 
models to bear in a phase- averaged formulation implies ail equivalence between 
ensemble- and spatial- or temporal- averages (ergodicity) that has been demon- 
strated neither experimentally nor analytically. Use of conventional models also 
demands a reasonable degree of commonality in turbulence structure between the 
benchmark flow' and the engine. While universality has been argued in the limit 
of fully- developed high-Reynolds-number broad- inertial-range turbulence, it is du- 
bious for the low Reynolds numbers (Section 4) and complex three-dimensional 
transient flows that characterize the engine. 

The same moderate Reynolds numbers that make IC engine flow problematic 
for RANS render it an appealing candidate for LES. It has been estimated that 
grid-independent (to a 10%-20% level) RANS computations of in-cylinder flow' and 
combust ion require at least 100 J mesh points using second-order or higher numerical 
methods. This corresponds to sub-millimeter mesh spacings in a typical automotive 
IC engine, and is not far beyond current practice of 250,000 to 500,000 nodes. This 
mesh density should suffice to capture large- and intermediate-scale flow structures. 
Thus for IC engines. LES mesh requirements are expected to be comparable to 
those of RANS. 

LES also promises more direct access to physical processes. Cvcle-to-cycle vari- 
ability in flow and combustion is one phenomenon that has proven elusive to ensemble- 
mean modeling and analysis. The result has been a number of largely ad hoc 
attempts to distinguish among ‘mean,’ 'turbulence/ and 'cyclic variability’ compo- 
nents of the flow (El Tahry Sc Haworth 1992). This distinction becomes moot in a 
spatially-filtered formulation. 

A drawback of LES is that substantially more computational effort may be re- 
quired compared to RANS. In the engine, for example, accumulating ensemble- mean 
statistics via LES requires computations through multiple engine cycles. Other is- 
sues include geometric complexity (moving piston and valves) and the relatively 
immature state of LES for modeling complex engineering flows. It is the purpose 
of the present research to advance LES on several of these fronts. First, we seek 
to establish the limitations and resolution requirements of a particular numerical 
methodology (Section 2). Second, w T e implement and evaluate the performance of 
state-of-the-art subgrid-scale turbulence models on unstructured deforming meshes 
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(Section 3). And third, we consider several other physical modeling issues that 
arise in engines, including treatment of solid walls and inflow/outflow boundaries 
(Section 4). 


2. Numerical methodology 

The high-order finite difference schemes and spectral methods that traditionally 
have been used for DNS and LES are not suited to complex geometric configurations 
with moving boundaries. On the other hand, the first-order time, second-order space 
discretizations typically employed for engineering RANS computations (Amsden et 
al. 1989, Haworth et al. 1990) are overly dissipative for LES. Here we require both 
that the numerical methodology be compatible with the intended application of 
in-cylinder flow and combustion, and that it be sufficiently accurate for meaningful 
large-eddy simulation. 

A novel discretization scheme called NO-UTOPIA ( NO de- Centered Unstructured 
TOpology, Parallel Implicit Advection) has been developed recently by O'Rourke 
& Sahota (1996a, 1996b). NO-UTOPIA is a variable explicit/implicit advection 
scheme that differences along characteristic directions to yield formal second-order 
accuracy in space and time, provided that the material-speed Courant number is 
less than unity. It has been implemented using node-centered variables and an 
edge-based data structure, allowing fully unstructured meshes. 

For the present study, the equations solved are those of conservation of mass, 
momentum, and enthalpy; the equation of state is that of an ideal gas with con- 
stant specific heats. Computations are compressible, although the Mach number 
is small for all cases considered here. To accommodate arbitrary mesh deforma- 
tion, advective terms in the governing equations are written using velocities relative 
to the moving grid. A pressure-based iterative solution procedure patterned after 
SIMPLE is used. 

The momentum equation for a control volume is obtained by integrating the 
Navier-Stokes equations over an arbitrary volume V with bounding surface 5 that 
is moving with velocity v. Density, Cartesian velocity components, and pressure 
are denoted by p, Uj, and p, respectively; the (constant) laminar viscosity is p L - 
Adopting Cartesian tensor notation with summation over repeated lower-case Ro- 
man indices, the momentum equation has the form, 


£ 

dt 


L 


puidr + 



— vj)dAj 


j pdA, + J T ef[ ],dA } + pg,dr . 


( 1 ) 


Here g , is a body force per unit mass. The effective stress r e r Jt is the sum of a 
viscous or laminar stress tl and the subgrid-scale stress tsgsji, 


7~eff ji 


Here S Jt is the rate-of-strain tensor, Sji = (§^- + §t|-)/ 2, and 6 Jt is Kronecker s 
delta. No supplementary turbulence model transport equations are carried for the 
subgrid-scale models considered here (Section 3). 
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Care is needed in the specification of initial conditions. For each flow, we select 
a reference velocity and length Uq and Lq (Section 4 ). Fluid properties 7 = c p /c v 
(ratio of specific heats) and R (specific gas constant) are set to nominal values for 
air. The Mach number based on U 0 is M 0 = U 0 /c Q where eg = 7 RT 0 is the square 
of the reference sound speed and To the reference temperature. Reference density 
and pressure are p 0 and po, respectively. Laminar viscosity pt is set to match the 
desired reference Reynolds number Re 0 . Remaining reference quantities and fluid 
properties are specified as: 

7 = 1.4 , R = 288.291 J/kg - K , M 0 = 0.1 , 

A, = 1.0 , po = (jMq ) -1 , = p 0 UoLoRe ^ . (3) 

Initial nodal velocities and pressures are prescribed, and nodal temperatures are set 
such that the initial entropy is uniform, T = 2o(p/po)* 7-1 ^ 7 . 


3. Subgrid-scale models 

The three models considered are of the Smagorinsky type. Here the influence of 
unresolved (subgrid-scale) motions on the resolved scales is treated as an additional 
viscosity, so that t$gs ji has a form identical to that of Jt , 

r SGS ji = 2psGS‘^>i ~ 2psGS &ji/ 3 . (4) 

The subgrid-scale viscosity psGS is taken to be proportional to a norm of the local 
rate-of-strain |S| and to a filter width A, 

^SGS = /JGA 2 |S| , with |5| = 2(S,jS,j) 1/2 . ( 5 ) 

The single model coefficient is C 3 . It is the specification of C 3 and A that distin- 
guishes the three models. 

3.1 Constant- coefficient Smagorinsky model 

The simplest model results from taking C 3 to be a constant and A to be equal 
to the mesh spacing. To accommodate non-uniform mesh spacing, A in Eq. ( 5 ) is 
specified as, 

A = V 1 / 3 , (6) 

where V is the volume associated with a computational element. 

Calibration with respect to benchmark turbulent flows has led modelers to adopt 
different values of C 3 . For homogeneous isotropic decaying turbulence, the value 
C 3 — 0.17 2 is found to result in a good match with the measurements of Comte- 
Bellot k Corrsin (1971) (Section 4.1). For planar channel flow, a value of C 3 = 0.1 2 
yields better agreement with measurements and DNS data (Section 4 . 3 ). The non- 
universality of C s motivates the need for a more general model. 
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8.2 Dynamic Smagorinsky model 

Germano et ai (1991) proposed an approach for evaluating subgrid-scale model 
coefficients from information contained in the resolved fields. Two filter widths are 
introduced, A and A, where A > A. Dependent variables filtered at scale A are 
denoted by the overbar notation () while the hat notation () denotes filtering at 
the larger scale. Formally, the first filter corresponds to an explicit filtering of the 
governing equations at the scale A. In practice, the first filter is implicit in the 
numerical method. That is, the quantity Ui(x,t ) denotes the computed velocity 
delivered by the numerical method at position x_ and time t . 

The second filter, or 'test filter,’ corresponds to a hypothetical second filtering at 
a larger scale. Thus Ui(x_,t) represents the LES-computed velocity field filtered at 

scale A. Stress tensors r Jt and Tj t represent the subgrid-scale stresses for the two 
filter widths, respectively: r Jt = pu 3 u t — pxijUi , and Tji = pu ) u l — pu 3 u x . Here the 
fluid density p can vary in time, but is nearly uniform in space (lo^Machjmmber). 
Filtering tj % through the second filter yields the quantity Tj t = pUjU % — pUjUt. Then 
subtracting r Jt from T Jt yields a second-order tensor Lji, which can be thought of 
as the stresses resulting from turbulent motions at scales intermediate between A 

and A: ^ ^ 

= Tji - Tji = pxijUi - pUjUi . (7) 

Equation (7) is the Germano identity. 

In the Smagorinsky model, subgrid-scale stresses at both scales are modeled con- 
sistently as, 

Tj i = 2pCsE 2 \S\Sji , and T Jt = 2 P cS . (8) 

Equation (8) is substituted into Eq. (7) to yield, 

Lji = 2pC a t 2 0 \%i - 2pC a A 2 \S\S ji . (9) 

Under the assumption that C s and A are constants with respect to the second (test) 
filtering operation, 

Lji = 2pC a A 2 ( (i/A) 2 |f|f J , - jSi^ ) = -2pC a A 2 Mj, . (10) 

Equation (10) defines the second-order tensor Mji which, like L Jt , is directly com- 
putable from the LES resolved velocity field. 

The quantity C a A 2 is chosen in a manner that minimizes the error in satisfying 

Eq. (10), t 3 i = L Jt + 2C a A 2 Mji. Here we follow Lilly (1992) in minimizing this 

2 

error in a least-squares sense with the constraint that C 3 A does not vary over 
homogeneous directions to yield, 


r A 2 — (TtjM^) 

“ 2 (M ki M k ,) ‘ 


( 11 ) 
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The angled brackets ( ) represent an average over homogeneous directions. 

The dynamic model offers several advantages compared to the const ant -coefficient 
model. Subgrid-scale viscosity increases locally in areas of low grid resolution in 
response to the high energy found between the two filter scales (large Lji). And, 
j^sgs decreases to zero in case all scales of motion are fully resolved locally ( L Jt — ► 0). 
A second advantage of the model as formulated here is that the filter width itself 

A need not be explicitly specified. It is the ratio of the filter widths A/ A that 
appears in M l} (Eq. 10). It is expected that the ratio of filter widths should be 
more uniform than the filter width itself on nonuniform deforming meshes. 

A fundamental limitation of the model as outlined here is that it requires at least 
one direction of statistical homogeneity. We therefore consider a third variant of 
the Smagorinsky model that addresses this shortcoming. 


3.3 Lagrangian dynamic Smagorinsky model 


Meneveau et al. (1996) proposed to accumulate the averages required in the 
dynamic model over flow pathlines rather than over directions of statistical homo- 
geneity. In this case, the error incurred by substituting the Smagorinsky model 
(Eq. 8) into the Germano identity (Eq. 7) is minimized along fluid-particle trajec- 
tories. The error to be minimized is the accumulated local squared error E along 
the pathline followed by the fluid particle that is located at position x at time t: 
E = /loo e nU where z(t') = x- J* u(z(t"),t")dt" is 
the trajectory followed by the fluid particle at times t* < t. The quantity W(t — t f ) 
is a weighting function that determines the relative importance of past events, and 
the error e t) is the difference between left- and right-hand sides of Eq. (10). As 

in the previous formulation (Section 3.2), it is assumed that C S A does not vary 

2 

strongly over the scale of the test filter. Then the value of C 3 A that minimizes the 
error E is, 


C, A 2 


1lm_ 

1mm 


( 12 ) 


where, 







LijMijUifMJWO-l'W , 




(13) 


An expedient choice of weighting function is one that decays exponentially back- 
wards in time, W(t — t') = T -1 exp [— (t — t')/T ], T being a memory or relaxation 
time scale. This choice allows the integral quantities li m and 1mm to be expressed 
as solutions to transport-relaxation equations. Meneveau et al. (1996) observed 
that high numerical accuracy is not needed in the solutions of these equations, and 
adopted the expedient of updating nodal values of Tim and 1mm by interpolating 
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from surrounding nodes at the upstream position: 

ntxM) 

where a = (At/T n )/(l + A t/T n ). Here superscript n + 1 denotes quantities eval- 
uated at the current time, superscript n quantities at the previous time, and A t is 
the computational time step. Trilinear interpolation is used to evaluate quantities 
at the upstream position x — ilAt from computed values at the surrounding nodes. 

We adopt the relaxation time T selected by Meneveau et al (1996), 

T = 9A[ (2A 2 ) 3 1 lm 1 mm T 1 * * * ' 8 , (15) 

where the value of the model coefficient is 0 = 1.5. This choice for T tends to reduce 
the memory time in regions of high strain (large and in regions of large 

nonlinear energy transfer (large Li 3 M XJ ). The memory time increases to reach back 
further in time along the particle’s trajectory in case L XJ M X3 remains negative over 
a persistent period. Negative L tJ M XJ might otherwise result in negative subgrid- 
scale viscosity, implying energy transfer from small to large scales and numerical 
instability. 

4. Flow configurations 

In a reciprocating engine, all flow velocities scale with the mean piston speed, 
which is proportional to the crankshaft rotational speed; length scales are inde- 
pendent of engine speed. Thus the mean-flow Reynolds number Reb (based on 
bore diameter and mean piston speed) and the turbulence Reynolds number Re\ 
(based on turbulence intensity and integral length scale) increase in proportion to 
engine speed. At 2,000 r/min, these are estimated to be Reb « 36,000 and Ret % 
1,000, respectively. In-cylinder turbulence, particularly at low engine speeds, is a 
low-to-moderate Reynolds number phenomenon. 

The number of turbulence ‘eddy-turnover’ times available for the decay of induction- 
generated turbulence in the engine is estimated to be greater than ten. Induction- 
generated turbulence has largely decayed by the time of ignition: it is the breakdown 
of the large-scale induction-generated flow structure that has the major influence on 
near-TDC turbulence and flame propagation. During compression and expansion, 
the in-cylinder flow is subjected to linear mean strains. The mean strain due to 
piston motion is slow compared to turbulence time scales, but persists for a large 
number of eddy-turnover times. These observations guide our choice of test cases 
for LES. 

4.1 Decay of homogeneous isotropic turbulence 

This canonical configuration is of relevance to the engine by virtue of the long 
times available for turbulence decay between intake-valve closure and ignition. 
Benchmark measurements were reported by Comte-Bellot & Corrsin (1971). There 


= max^ 0 , a[L t) M l j] rx+l {x ) + (1 - a)ll M {x - u n At) J , and 
= a[M tj M t j] n +\x) + (1 - a)ll fM {x - u n At) , 
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the temporal decay of homogeneous isotropic turbulence was approximated by grid- 
generated turbulence in a wind tunnel. 

Here we compare computed turbulence kinetic energy decay k = lfu x u l )j2 and 
three-dimensional energy spectra to the experimental data of Comte- Bellot &; Corrsin 
(1971). Turbulence was generated using a grid spacing of M = 5.08 cm in a uni- 
form mean flow of velocity U Q 0 = 10 m/s, yielding a Reynolds number of U^M/u 
= 34,000. Data were reported at three downstream stations: Ur^t/M = 42, 98, and 
171. 

Computations are done on triply-periodic uniform cubic meshes of length 27 r along 
each edge. Scalings are such that the edge of the computational box 2n corresponds 
to a physical length of 0.55 m, and the computational reference velocity Z7 0 = 1.0 
corresponds to the physical velocity Uoo = 10.0 m/s. Other scalings and parameters 
are summarized in Eq. (3). 

The initial velocity field is prescribed by a procedure similar to that used for 
incompressible spectral simulations. We begin with a superposition of Fourier modes 
having a prescribed energy spectrum but random phases; this is projected onto 
the divergence-free space. The resulting field represents the flow upstream of the 
first measurement station. It is advanced in time over several turbulence eddy 
turnover times to adjust to compressibility and to build phase coherence. The 
process is repeated with different initial fields until a satisfactory match is obtained 
between the computed and measured energy spectrum at the first measurement 
station U^tjM = 42. 

Comparisons between model and measurement are made on the basis of filtered 
quantities. Energy spectra are filtered based on our limited understanding of the 
nature of the filtering implied by the numerical method. We assume that () corre- 
sponds to a top-hat filter in physical space at the mesh spacing with trapezoidal-rule 
integration. 

Initial computations are for a mesh of 32 3 nodes. This has been the traditional 
starting point for new numerical methodologies in LES, but is marginal for resolving 
the physics of the flow. At the initial measurement station U^t/M = 42, the 
computational box edge corresponds to between ten and twenty integral length 
scales of the turbulence: fewer than three mesh points span an integral scale. By 
the final measurement station Uoot/M = 171, the turbulence integral scale has 
roughly doubled. The computational time step is prescribed such that material 
Courant numbers are less than unity. 

4-2 Linear strain of initially isotropic turbulence 

Homogeneous strain of initially homogeneous isotropic turbulence has been a 
second canonical configuration for analysis, modeling, and experiment. Here we 
consider the linear expansion and the linear compression for their particular rele- 
vance to the IC engine. Results are compared to rapid-distortion theory (RDT), 
a linearized theory of turbulence that is appropriate in the limit where the mean 
strain rate is large compared to an inverse turbulence eddy turnover time (Kassinos 
Sz Reynolds 1994). Although IC engines appear to be far from the RDT limit, this 
nonetheless provides a sound basis for model evaluation. This configuration also 
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exercises the code’s mesh deformation capability. 

The behavior that we seek to capture is the distribution of energy among the three 
normal-stress components. We monitor the evolution of the normalized anisotropy 
tensor b l} as a function of the total strain C*, 

hj = {u t u } )/(u,u,) - Sij/Z , C* = exp[ J \S*\dt ] . (16) 

Here S* is the dominant eigenvalue of the modified rate-of-strain tensor S* r where 
S*j = S tJ — SuSij! 3. In the absence of mean rotation, the evolution of b tJ (C*) 
in the RDT limit for a non-solenoidal mean strain Sij is the same as the evolu- 
tion of bij(C*) for the corresponding divergence-free rate-of-strain S*j (Kassinos &: 
Reynolds 1994). 

The linear expansion is the superposition of spherical (isotropic) expansion and 
irrotational axisymmetric contraction. Experiments (see Kassinos & Reynolds 1994 
for references) show that the anisotropy b tJ depends weakly on the magnitude of the 
mean rate of strain. Thus even for slow linear expansions, the evolution of hiflC*) 
is expected to be similar to that predicted by RDT. The linear compression is the 
superposition of spherical compression and irrotational axisymmetric expansion. In 
this case, experiments reveal that stronger anisotropy develops at slower mean rates 
of strain. 

Initial meshes and velocity fields are the same as those for decaying turbulence 
simulations (Section 4.1). The mean strain rate is imposed by deforming the domain 
in a manner that maintains a constant rate-of-strain 533 along the £ 3 coordinate 
direction. The mesh deformation rate varies linearly from zero at x 3 = 0 to SwL$(t) 
at £3 = Z 3 (f), yielding exponential expansion or contraction of the mesh with time, 

1.3(0 = £3(0) * exp[S 33 t]. 

4.8 Planar channel flow 

The planar channel flow adds the complexity of walls and a single statistically 
nonhomogeneous direction. Computations are performed on a rectangular prism of 
dimension L\ (streamwise) by L 2 (normal to the wall) by L 3 . Relevant dimensionless 
parameters are Reynolds numbers based on the wall friction velocity t/ r , and on the 
bulk velocity: Re r = u r 8/v, and Res = Ub^/v where Ub = J 0 2 (ui (£ 2 ))d£ 2 /T 2 , 
and 8 is the channel half-width. Here angled brackets ( ) denote averages over 
planes parallel to walls. 

Results are computed for a low Reynolds number of Re r = 180 (Res ~ 2,800). 
The computational domain is Ait 8 by 28 by 47r6/3. The initial mesh of 33 x 65 
x 33 nodes is comparable to that adopted by earlier researchers for this Reynolds 
number, although higher-order numerical methods have been used in most previous 
work (e.g., Piomelli 1993). Mesh spacing is uniform in x\ and £3 and follows a 
tanh distribution in £ 2 . Grid spacing varies from a minimum of A = 0.87 at the 
wall to a maximum of A = 11.7 at the channel centerline, where the + notation 
denotes standard wall-units scaling (y+ = yu T /v). Computations are periodic in 
and £3, with no-slip boundaries at £ 2 = 0 and £ 2 = i 2 . 
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FIGURE 1. Decay of filtered turbulence kinetic energy for homogeneous isotropic 
decaying turbulence. Filter corresponds to a top-hat in physical space on a 32* 
uniform mesh with trapezoidal-rule integration. Symbols are measurements of 

Comte-Bellot. & Corrsin (1971): o . Lines are computations: no subgrid-scale 

model; constant-coefficient Smagorinsky model with C s = 0.17 2 ; dy- 

namic Smagorinsky model. 

Periodicity in the streamwise direction is maintained by imposing a streamwise 
body force g x (Eq. 1) consistent with the desired Re r . A global force balance in the 
x i direction yields, 

9\ = 2L~ i {isRe r /fi) 2 . (17) 

Fluid viscosity is set to maintain the desired bulk Reynolds number Rcq. Velocity 
and length scales are t o = U b and L 0 = <5; remaining parameters are set accord- 
ing to Eq. (3). The flow is allowed to develop for about 20 flow-through times 
^dev t' b ! L\ zz 20 (t<i ev Ur/h ~ 15). Profiles of mean velocity. Reynolds stresses, and 
other statistics as functions of x <2 then are accumulated by averaging over planes 
parallel to walls and over time using averaging times f avg of at least t^ vg u T /S = 5. 
Computed results are compared to DNS results of Kim et al (1987) at the same 
Re r . 

4-4 Axisymmetric piston- cylinder assembly 

The target configuration for establishing the feasibility of in-cylinder LES is the 
simplified piston-cylinder assembly of Morse et al. (1978). There flow enters a 
pancake (flat head and piston) chamber through a central pipe of inside diameter 
18.75 mm and length 1.8 in. The piston is driven in simple harmonic motion at a 
speed of 200 r/min through a 60-mrn stroke; there is no compression. Bore diame- 
ter is 75- mm bore, and top- dead-center clearance height is 30-mm. Laser-Doppler 
anemometry has been used to obtain ensemble- mean (phase- aver aged) radial pro- 
files of mean and rms axial velocity at 10-mm axial increments starting at the 
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kLj 2n 

FIGURE 2 . Evolution of filtered three-dimensional energy spectra for homogeneous 
isotropic decaying turbulence. The filter is a top-hat of width 27r/16 in physical 
space with trapezoidal-rule integration. Symbols are experimental measurements of 
Comte-Bellot &: Corrsin (1971): • Uoot/M — 42; * Uoot/M — 98; + Uoot/M = 171. 
Lines are computations for the constant-coefficient Smagorinsky model with C 3 = 
0.17 2 on 32 3 meshes: Uoot/M = 42; Uoot/M — 98; Uoot/M = 171. 

head for crank positions of 36°, 90°, 144°, and 270° after top-dead-center. This 
flow can be thought of as an extension of the classic statistically stationary sudden 
expansion/contraction to a statistically periodic case. 

Several pieces of information are sought from these computations. First, we 
can evaluate the performance of subgrid-scale turbulence models in a configuration 
approaching that of an engine on a deforming unstructured mesh. Second, we will 
build experience with explicit phase- (ensemble-) averaging compared to spatial 
filtering and traditional RANS modeling. Third, we can establish mesh resolution 
requirements, particularly in the vicinity of walls. This includes a determination of 
the need for explicit wall models beyond that provided by the subgrid-scale model. 
And fourth, we will explore the nature of inflow forcing required to generate realistic 
in- cylinder flow variability. 

5. Results 

All displayed results represent the resolved motion delivered by the numerical 
method in combination with the specified subgrid-scale model. These are the ()- 
filtered quantities as defined in Section 3. No attempt has been made to add explicit 
subgrid-scale contributions to the stresses. 

5.1 Decay of homogeneous isotropic turbulence 

The effect of filtering on the fraction of resolved turbulence kinetic energy in the 
experiments of Comte-Bellot &: Corrsin (1971) has been computed. For the filtering 
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FIGURE 3 . Evolution of the normalized anisotropy tensor b t j as a function of 
total strain C* for the linear expansion. Open symbols are RDT results (Kassinos 
& Reynolds 1994 ): o RDT, 6n; A RDT, 622 \ D RDT, 633. Lines are computations 
for the constant-coefficient Smagorinsky model with C s = 0 . 17 2 on 32 3 meshes: 

~ S33 * ~ 4 , 611; 5*33 * k/e % 4 , 622; ^ S33 • k/e % 4 , 633; 533 • 

k/e % 8, 6n; 533 • k/e & 8, 622; 4 ++ 1 I S33 * k/e % 8, 633. (Results for the 

lower rate-of-strain S 33 are indistinguishable from those at the higher S33.) 

assumed to be closest to our numerical method (top-hat filter with trapezoidal- 
rule integration) only about 45 % of the energy is resolved at the first measurement 
station on the 32 3 mesh. 

The decay of filtered turbulence kinetic energy versus time for the 32 3 mesh is dis- 
played in Figure 1. With no subgrid-scale model, there already is substantial decay 
resulting mainly from numerical dissipation. Constant-coefficient Smagorinsky adds 
sufficient additional dissipation to yield good agreement with measurements, using 
the standard value of the model coefficient (C 3 = 0 . 17 2 ). The dynamic Smagorinsky 
model yields similar results, returning a value of C s % 0 . 16 2 , close to the standard 
value. 

Filtered three-dimensional energy spectra are plotted in Figure 2 . There is a 
pile-up of energy at wave numbers just beyond the peak of the spectrum in the 
computations. Thus while we are able to match the energy decay rate on this coarse 
mesh, the dynamics of the system are not fully captured. This is not surprising in 
a computation where less than half of the energy is resolved. 

5.2 Linear strain of initially isotropic turbulence 

Evolution of the normalized anisotropy tensor as a function of total strain is given 
in Figs. 3 and 4 . Results are presented for two different values of S33 • k/e to show 
the influence of mean rate-of-strain. All numerical results are for a 32 3 mesh using 
the constant-coefficient Smagorinsky model (C 9 = 0 . 17 2 ). RDT data are shown for 
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Figure 4. Evolution of the normalized anisotropy tensor b t j as a function of total 
strain C * for the linear compression. Open symbols are RDT results (Kassinos & 
Reynolds 1994): o RDT, b u ; * RDT, b 2 2 \ o RDT, 633- Lines are computations 
for the constant-coefficient Smagorinsky model with C, = 0.17 2 on 32 3 meshes: 

S33 ~ ~ 4 , &11; ~ - 4 , 622; S33/ fc/e « -4 , 6 33 ; 

S33 ■ k/e % - 8 , 6 n ; S33 ■ k/e » - 8 , b 22 ; -H-H+ S33 • k/e « - 8 , 6 33 . 


comparison. 

For the linear expansion, computations are in good quantitative agreement with 
RDT and are insensitive to the applied mean rate-of-strain (Fig. 3). This is consis- 
tent with experimental trends reviewed by Kassinos & Reynolds (1994). 

Results for the linear compression warrant further discussion (Fig. 4). In this 
case computed results are closer to RDT for the slower mean rate-of-strain, and the 
degree of anisotropy increases with increasing mean rate-of-strain. This is contrary 
to experimental trends, which show increasing anisotropy at slower rates of strain 
(Kassinos &; Reynolds 1994). 

The Reynolds-averaged turbulence stress transport equation for homogeneous 
turbulence subjected to a uniform mean strain rate is derived by standard proce- 
dures, 


d(p{u' k u\)) 


d(u,) 


dt 


= i> 


j—t \ &{ u k) 


dxi 


+ T ki + T kl - tkt • 


( 18 ) 


The prime notation emphasizes that there is non- zero mean flow, U' = u l - (u t ). 
Here the first two terms on the right-hand side represent the rate of production, 
and T* kl are the ‘rapid 7 and ‘slow 1 pressure-rates-of-strain, respectively, and 
is the viscous dissipation. In the limit of rapid distortion, e k i and T kl are negligible. 

For the present linear mean strains, d{u x )jdxj = 533^,3^3. Thus all turbulence 
production goes directly into (u l 3 2 ) and is redistributed to the other two components 
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FIGURE 5. Streamwise mean velocity profiles normalized by the bulk velocity 
Ub for the planar channel flow at Re r = 180. Symbols ( * ) are the DNS data of 

Kim et al (1987). Lines are computations on the 33x65x33 mesh: constant- 

coefficient Smagorinsky, C 9 = 0.17 2 ; dynamic Smagorinsky; Lagrangian 

dynamic Smagorinsky. 

(u\ 2 ) and {u' 2 2 ) via the pressure-rate-of- strain terms. For the low-resolution LES 
computations of linear compression, the effective rapid pressure-rate-of-strain model 
does not redistribute sufficient energy from the ‘direct’ production component to 
the other two. Moreover, the effective slow pressure-rate-of-strain model responds 
incorrectly to a decrease in the mean rate-of-strain. 

5.S Planar channel flow 

Mean velocity profiles from the dynamic Smagorinsky and Lagrangian dynamic 
Smagorinsky models are very similar to one another, and show better agreement 
with DNS than the constant- coefficient model (Fig. 5). All three models deviate 
from DNS in the logarithmic region (y+ > 10). Ratios of centerline mean velocity 
to bulk velocity (ui(y = S ))/Ub are 1.22 for constant-coefficient Smagorinsky, 1.15 
for dynamic Smagorinsky, 1.15 for Lagrangian dynamic Smagorinsky, and 1.16 for 
the DNS of Kim et al. (1987). 

Both dynamic models effectively ‘turn down’ the subgrid-scale viscosity in the 
vicinity of the wall. The mesh spacing A decreases close to the wall, as does 
the model coefficient C s . The latter behavior is demonstrated in Fig. 6. There 
computed profiles of CJ extracted from the dynamic model and the Lagrangian 
dynamic model are shown. For the former model, the standard value of 0.1 is 
recovered in the center of the flow, with a rapid drop-off to zero at the walls. The 
Lagrangian dynamic model behaves similarly out to a distance of about y + « 40, 
but levels off to a lower value of C\^ 0.06 at the centerline. 

Computed Reynolds-stress profiles from the Lagrangian dynamic model are given 
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FIGURE 6. Computed profiles of C 9 adjacent to the lower wall for the planar 
channel flow at Re r = 180. Symbols (•,*,+) are planar* averaged profiles from the 

dynamic Smagorinsky model at three instants of time. Lines ( , , ) 

are planar- averaged profiles from the Lagrangian dynamic Smagorinsky model at 
three instants of time. 


in Figs. 7 and 8. Results from the dynamic model are similar, while the constant- 
coefficient model yields somewhat poorer profiles (not shown). This is consistent 
with our findings from the mean velocity profiles of Fig. 5. Normal stress com- 
ponents display qualitatively correct behavior (Fig. 7), but there are significant 
quantitative departures from the DNS data. In particular, on this coarse mesh, all 
models tend to leave too much energy in the direct production component (u\ 2 ) 
at the expense of (u^ 2 ) and (u^ 2 ). The value of the peak shear stress is computed 
reasonably well, although the LES profile is shifted outward from the wall compared 
to the DNS data (Fig. 8). These findings suggests that the present mesh resolution 
is marginal for computing second-order statistics, especially in the log-law region. 

5.4 Axisymmetric piston- cylinder assembly 

Computations are in progress at the time of this writing. Quantitative compar- 
isons with measurements of Morse et al. (1978) are forthcoming. 

6. Discussion and conclusions 

This research has explored a candidate numerical methodology and subgrid-scale 
stress models for LES of flow in reciprocating IC engines. The present results 
have been obtained using coarse meshes that are representative of minimal mesh 
requirements for spectral LES. Generally reasonable evolution of first and second 
moments has been found nevertheless. This is an encouraging finding, given the 
low formal accuracy of the numerics. Based on these early results, it is anticipated 
that acceptable accuracy can be obtained using practical mesh densities. 
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FIGURE 7. Turbulence intensities normalized by the wall friction velocity u T for the 
planar channel flow at Rt r = 180. Symbols are the DNS data of Kim et ai (1987): 
• streamwise (xi) component; * wall-normal (# 2 ) component; + cross-stream ( 2 : 3 ) 
component. Lines are computations using the Lagrangian dynamic Smagorinsky 

model (resolved portion): streamwise (xi) component; wall-normal 

(£ 2 ) component; cross- stream (£ 3 ) component. 

Specific deficiencies have been attributed to inadequate spatial resolution. These 
include the energy spectrum decay for isotropic turbulence and insufficient energy 
transfer from the 'direct 1 production component for linear compression and planar 
channel flow. The two dynamic models have demonstrated an advantage compared 
to the constant-coefficient model in the planar channel flow. No specific deficiencies 
of the dynamic subgrid-scale models have been identified. In some cases, model 
results are not much different than those obtained in the absence of any explicit 
subgrid-scale model. This is consistent with earlier LES work for coarse meshes 
and low-order numerical methods. It remains to establish that these deficiencies 
can be overcome through mesh refinement, and to quantify resolution requirements 
for a specified level of fidelity to experiment or to benchmark computations. Short 
of explicitly filtering the governing equations at a scale much larger than the mesh 
spacing, it will remain difficult to isolate numerical inaccuracy from subgrid-scale 
model performance in LES. 

Beyond spatial resolution, the most pressing outstanding issue is the lack of an- 
alytic characterization of the filtering implied by non-spectral numerical methods: 
what is U, ? While it is straightforward to analyze and implement a variety of filters 
in spectral methods (e.g., spectral cutoff, spatial top-hat, spatial/spectral Gaus- 
sian), there has been little analysis to guide the implementation of filters implicit in 
finite-difference, finite- volume, or finite-element schemes on unstructured meshes. 
Our experience with the initial spectrum for decaying turbulence shows that the 
present discretization scheme affects all wavenumbers to some extent. The same has 
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FIGURE 8. Turbulence shear stress normalized by the square of the wall friction 
velocity for the planar channel flow at Re r — 180. Symbols (o ) are the DNS data of 

Kim et al . (1987). Lines ( ) are computations using the Lagrangian dynamic 

Smagorinsky model (resolved portion). 

been found for other non-spectral methods that are being explored for unstructured 
LES (Jansen 1995). 

Other outstanding issues for in-cylinder LES include wall treatment and inflow 
boundary conditions. Piomelli (1993) has shown that accurate LES results can 
be obtained using the dynamic model at high Reynolds numbers without further 
explicit wall modeling. The challenge at inflow boundaries is to establish the nature 
of forcing needed to yield in-cylinder velocity statistics representative of measured 
‘cyclic variability. 1 A final determination of suitability awaits the results of finer- 
mesh simulations for the three canonical configurations, and multiple-cycle results 
for the axisymmetric piston-cylinder assembly. 
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